Modulational Instability and Complex Dynamics of Confined Matter- Wave Solitons 
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We study the formation of bright solitons in a Bose-Einstein condensate of 7 Li atoms induced by a 
sudden change in the sign of the scattering length from positive to negative, as reported in a recent 
experiment (Nature 417, 150 (2002)). The numerical simulations are performed by using the 3D 
Gross-Pitaevskii equation (GPE) with a dissipative three-body term. We show that a number of 
bright solitons is produced and this can be interpreted in terms of the modulational instability of the 
time-dependent macroscopic wave function of the Bose condensate. In particular, we derive a simple 
formula for the number of solitons that is in good agreement with the numerical results of 3D GPE. 
By investigating the long time evolution of the soliton train solving the ID GPE with three-body 
dissipation we find that adjacent solitons repel each other due to their phase difference. In addition, 
we find that during the motion of the soliton train in an axial harmonic potential the number of 
solitonic peaks changes in time and the density of individual peaks shows an intermittent behavior. 
Such a complex dynamics explains the "missing solitons" frequently found in the experiment. 
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The recent experimental observation of dark [1] and 
bright [2,3] solitons in Bose-Einstein condensates has re- 
newed the interest in the intriguing dynamical phenom- 
ena of nonlinear matter waves [4]. In the experiment 
of Strecker et al. [3] soliton trains have been formed by 
making a stable condensate of 7 Li atoms with a large 
positive scattering length a s using a Feshbach resonance 
and then switching a s to a negative value. The forma- 
tion of the soliton train has been interpreted as due to 
quantum mechanical phase fluctuations of the Bosonic 
field operator [5]. By imposing a suitable space depen- 
dent pattern in the initial phase of the Bosc condensate 
and then using the ID time-dependent Gross-Pitaevskii 
equation (GPE), Al Khawaja et al. [5] have reproduced 
the formation of the soliton train. 

In this paper we show that the formation and sub- 
sequent evolution of a soliton train can be adequately 
investigated by using the classical (mean-field) time- 
dependent 3D GPE with a dissipative term which takes 
into account the three-body recombination process that 
is crucial during the collapse of the condensate [6]. The 
multi-soliton configuration is obtained without imprint- 
ing the initial wave function with a fluctuating phase. 
We show that the soliton train is produced by the mod- 
ulational instability (MI) of the evolving classical phase 
of the Bose condensate (see also [7]). MI is a nonlin- 
ear wave phenomenon in which an exponential growth of 
small perturbations results from the interplay between 
nonlinearity and anomalous dispersion. MI has been pre- 
viously studied for waves in fluids [8], in plasma physics 
[9], in nonlinear optics [10], and also in the context of 
the superfluid-insulator transition of a Bose-Einstein con- 
densate trapped in a periodic potential [11]. Here we find 
that the number of bright solitons induced by MI is given 
by a simple analytical formula which reproduces our nu- 
merical simulations. By investigating the long time evo- 
lution of the soliton train under the action of a harmonic 



potential of frequency lo z in the travelling axial direc- 
tion we find that the center of mass motion is periodic 
with frequency oj z but the density of each solitonic peak 
strongly changes in time. 

At zero temperature the macroscopic wave func- 
tion i/j(r, t) of a Bose-Einstein condensate made of 7 Li 
atoms can be modeled by the following dissipative time- 
dependent Gross-Pitaevskii equation 



ih^- + ^V 2 -U- gN\if + ijN 2 ]^ 4 
at 2m 



^ = 0, (1) 



where g = Airft 2 a s /m, a s is the s-wave scattering length, 
m the atomic mass, N is the number of condensed atoms 
and 7 is the strength of the dissipative three-body term 
[6]. At t = the condensate wave function is normalized 
to one. Following the experiment of Strecker et al. [3], 
the external potential U(r) is given by 

U(r) = | [ul (x 2 + y 2 ) + x 2 c 2 z 2 ] + V L {z) (2) 

where Vt, (z) is the box optical potential that initially con- 
fines the condensate in the longitudinal direction. The 
harmonic confinement is anisotropic with lo±_ = 2n x 800 
Hz and uj z = 2tt x 4 Hz. \ is a parameter that modifies 
the intensity of confinement in the axial direction. 

To investigate the formation of the train of bright 
solitons, we first calculate the ground-state wave func- 
tion of the condensate with a positive scattering length 
a s = lOOas, where as = 0.53A is the Bohr radius, by 
using Eq. (1) with imaginary time and 7 = 0. The 
numerical code implements in cylindric symmetry (z, r) 
a finite-difference splitting method (spatial grid with 
400 x 100 points) with a predictor-corrector algorithm 
to treat the nonlinear term [12]. We choose a conden- 
sate with longitudinal width L = 284.4 /im and iV = 10 4 
atoms. Then we use the ground-state wave function as 
initial condition for the time-evolution of Eq. (1) with 
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a s — — Sets, Vl{z) = and \ = 0. Following [6] we 
choose 7 = 1.77 x 1CP 11 {huj z )/a Q z . 7 is important during 
the collapse, because is fixes the critical density at which 
the compression ceases. 




FIG. 1. Axial density profile p{z) of the Bose-Einstein con- 
densate made of 10 4 7 Li atoms obtained by solving dissipative 
3D GPE. For t < the scattering length is a 3 = lOOas, while 
for t > we set a s = — 3a b with as the Bohr radius and 
Vl = 0. External harmonic potential given by Eq. (2) with 
X = 0. Length z in units a z — (h/mujz) 1 ^ 2 , density p in units 
l/a z and time t in units lo^ 1 . 



In Figure 1 we plot the probability density in the longi- 
tudinal direction p(z) = J dxdy\ip(x, y, z)\ of the evolv- 
ing wave-function. The initially homogeneous conden- 
sate shows the formation of 5 peaks. Note that ini- 
tially the phase of the condensate is set equal to zero. 
After the formation, the peaks start to separate each 
other showing a repulsive force between them. Even- 
tually each peak evolves with small shape oscillations 
but without appreciable dispersion. This means that 
strictly speaking the train is not made of shape-invariant 
waves. We will use the word "solitons" to indicate the 
peaks of the train because they remain spatially local- 
ized during the time evolution. The formation of these 
solitons can be explained as due to the modulational in- 
stability (MI) of the time-dependent wave function of the 
Bose condensate, driven by imaginary Bogoliubov exci- 
tations [5]. The Bogoliubov elementary excitations e k 
of the static Bose condensate $(r) can be found from 
Eq. (1) by looking for solutions of the form vp(r 7 t) — 
e-^/R [$( r ) +u k (r)e- l£kt / h + v* k (r)e l£kt / h ] , and keep- 
ing terms linear in the complex functions u(r) and v(r). 
Neglecting the dissipative term 7, in the quasi-lD limit 

[13] one finds e k = ^h 2 k 2 /(2m) (h 2 k 2 /(2m) + 2g 1D n), 

where gm = g/{2^a 2 \_) with aj_ = \fh] (nuJ±), n = N/L 
is the linear density and L is the length of the conden- 
sate. By suddenly changing the scattering length a s 
to a negative value, the excitations frequencies corre- 
sponding to k < k c — ^/167r|gi£>|n become imaginary 



and, as a result, small perturbations grow exponentially 
in time. It is easy to find that the maximum rate of 
growth is at fc = k c / \/2. The wavelength of this mode 
is A = 27r/fc and the ratio L/X gives an estimate of 
the number N s of bright solitons which are generated: 
N s = y/N\^\L/(ira±). 
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FIG. 2. Axial density profile p(z) of the Bose-Einstein con- 
densate made of 10 4 7 Li atoms at t — 4.6 obtained by solving 
dissipative 3D GPE for three values of the final scattering 
length a a . External harmonic potential given by Eq. (2) with 
X = 0. Scattering length a s in units of the Bohr radius as. 
Units as in Figure 1. 



As predicted by this formula, Figure 2 shows that 
the number of peaks grows with the scattering length. 
The predicted number N s of solitons is in very good 
agreement with the numerical results: N s = 2.60 for 
a s = —las, N s = 4.51 for a s = — 3clb and N s = 5.82 
for a s = —5as- The period tq = h/Im(ik ) associated 
to the most unstable mode is given by r = mL 2 / (ttHN 2 ) , 
in rough agreement with the numerically estimated for- 
mation time of the strongest fluctuations in the number 
of peaks belonging to the soliton train. 

In the experiment of Strecher et al. [3] there are ini- 
tially about 10 5 atoms. In this case our numerical (see 
below) and analytical results predict the formation of 
about N s — 15 solitons, that is twice the number exper- 
imentally detected. We have verified that the number 
of solitons does not depend on the dissipation constant 
7 of Eq. (1) but increasing 7 the densities of solitons 
are reduced and their widths are increased. Note that in 
this system recombination processes have not been inves- 
tigated experimentally and these phenomena can affect 
the imaging process of soliton train densities. 

We have also verified that the number of solitons in- 
creases by increasing At, the delay time between the re- 
moval of endcaps and the change of scattering length (in 
figure 1 and 2 it is At = 0). For instance, with the data 
and units of Figure 2, A^ s = 6 for At = 0.3 and N s = 7 for 
At = 0.6. That is a simple consequence of the enlarge- 
ment of the axial width of the Bose-Einstein condensate, 
in full agreement with the experimental results [3]. 
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FIG. 3. Ratio Nr/N as a function of time for two neigh- 
boring solitons with an initial phase difference <f) and spatial 
separation 2zo = 6. N = Nr + Nl is the sum of the number 
Nr of atoms in the right soliton and the number Nl of atoms 
in the left soliton. N\a s \/a± = 1/%/TO. Results obtained 
by solving conservative NPSE. External harmonic potential 
given by Eq. (2) with ujj_/uj z = 10 and \ — 0. Units as in 
Figure 1. 

In order to clarify the role of phases in the interaction 
between bright solitons we consider a simpler model: two 
solitons without dissipation. It has been found that the 
effective interaction between two bosonic matter waves 
depends on their phase difference <f> , being proportional 
to a s cos (<?!>) [5,14]. The axial wave function f(z) of 
a bright soliton under transverse harmonic confinement 
can be analytically determined [14] by using the non- 
polynomial Schrodinger equation (NPSE), an effective 
ID equation we have recently derived from the 3D GPE 
[15]. In [14,15] we have shown that ID GPE reproduces 
bright solitons of 3D GPE only when the Bose conden- 
sate is strongly cigar-shaped. Instead, NPSE always 
reproduces 3D GPE solitons with great accuracy. As 
initial condition for the numerical solution of the time- 
dependent NPSE (spatial grid of 10 4 points) we choose 
ip(z) = [f(z - zo) + f(z + z )e 1 *] /V2, where cj> is the 
phase difference of the two neighboring solitons centered 
in — zq and zq. 

In Figure 3 we plot the time evolution of Nr/N, where 
Nr is the number of atoms in the right soliton, for some 
values of <fi with N\a s \/a± = 1/VlO, ^±/^z = 10 and 
X = Q- For = the two solitons are attractive and 
eventually form a static peak which radiates small waves. 
For <f> — 7r/4 the centers of the two solitons do not change 
with time but they exchange atoms as in a Josephson 
junction. Obviously the amount of atom exchange will 
depend on the details of the two solitons at initial time 
(widths and separation distance). We find such a com- 
plex exchange of atoms for < cj) < 7r/2, while for 
<f> = tt/2 the two solitons are weakly repulsive, their 
shapes change with time and eventually two solitons of 
different density appear because the Josephson exchange 
cannot continue when the two solitons separate due to 
their repulsive interaction. Note that with <fi = —tt/2 the 
densities of the two peaks are inverted. For larger values 
of <j) the dynamics is quite similar to the (f> — tt/2 case, 
with two solitons of different shapes getting away from 



each other. The fraction of exchanged atoms between 
the two solitons goes to zero as 4> approaches tt due to 
the increase of repulsive interaction and the decrease of 
interaction time. Exactly at <j) = tt parity symmetry of 
the problem inhibit atom exchange. 
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FIG. 4. Moving soliton train made of Bose condensed 7 Li 
atoms obtained by solving the dissipative ID GPE. Initially 
there are N = 10 5 atoms. Three frames: t = 0.8 (top), t — 1.6 
(middle), and t = 3.8 (bottom). External harmonic potential 
given by Eq. (2) with \ = 1- Red color corresponds to highest 
density. Units as in Figure 1. 

In the soliton train shown in Figure 1 and Figure 2 
the phase difference of adjacent solitons is not imposed 
a priori at t — 0, like in Ref. [5], but is self consistently 
determined by the GPE. As a consequence, its value is 
not exactly equal to tt and it changes with time. It is 
interesting to study the long time evolution of the soli- 
ton train under axial harmonic confinement as done in 
the experiment of Strecker et al. using 10 5 7 Li atoms [3]. 
Instead of 3D GPE in this numerical simulation we use 
the ID GPE (spatial grid of 10 4 points) with rescaled ef- 
fective interaction and dissipative term [15], which allows 
to extend our calculations to long times. As explained 
in [15], ID GPE is obtained from 3D GPE by imposing 
a Gaussian wave function of width a± in the transverse 
direction. Note that we have verified that the ID GPE 
results are in good agreement with 3D GPE ones for the 
set of parameters of Streker et al. experiment [3] . 

In Figure 4 we show three frames of the traveling train 
as color contour-plots of the density. In each pannel the 
orizontal axis is the z coordinate and the vertical axis 
is the x coordinate. The center of mass of the soliton 
train oscillates with the frequency of the harmonic con- 
finement. Moreover, the solitons spread out in the middle 
and bunch at the turning points in very good agreement 
with the experimental results under similar conditions 
[3]- ^ ' 

It is interesting to observe that during the time evolu- 
tion the densities of bright solitons oscillate in an irregu- 
lar way and, at certain instants, a few solitons practically 
disappear and reappear after a while. Such intermittent 
behavior is related to the complex dynamics of the bright 
soliton phases due to the presence of axial harmonic con- 
finement. As confirmed in Figure 5, in absence of axial 
confinement (x = 0), after a transient the peaks, separat- 
ing each other, become true solitons. In particular, the 
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number N p of peaks in the train gets constant and equal 
to the analytically estimated number N s of bright soli- 
tons. Instead, in presence of axial harmonic confinement 
(x = 1), N p changes in time. Obviously, in this case the 
peaks never become true solitons. 
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FIG. 5. Number N p of peaks in the train of Bose con- 
densed 7 Li atoms as a function of time t. N is the initial 
number of atoms. Simulation performed by using the dissipa- 
tive ID GPE (dotted line). The solid line is obtained from a 
smoothed finite-resolution density profile of the train. Exter- 
nal harmonic potential given by Eq. (2). Units as in Figure 
1. 



solitons is overcome by the potential energy of solitons 
and the solitons cross in the middle of the trap. We have 
verified that to get spreading solitons in the middle of the 
trap with initially 10 4 atoms it is necessary to strongly 
reduce the axial frequency (x = 1/4). 

In conclusion, we have successfully explained and nu- 
merically simulated the dynamical process of soliton train 
formation induced by modulational instability in the 
framework of the time-dependent Gross-Pitaevskii equa- 
tion with a three-body dissipative term. Contrary to 
the claim of Al Khawaja et al. [5], it is not necessary 
to include quantum phase fluctuations to trigger the for- 
mation of the soliton train. We have also investigated 
the soliton-soliton interaction and found a novel phe- 
nomenon: the intermittent dynamics of individual peaks 
during the time evolution of the soliton train in an axial 
harmonic potential. Signatures of this phenomenon can 
be extracted from the data of the experiment of Strecker 
et al. [3]. Because of the intimate connection between 
atom optics with Bose-Einstein condensates and light op- 
tics we believe that the intermittent dynamics in soliton 
trains may be also observed with optical solitons in fibres. 

L.S. thanks R.G. Hulet for many useful e-discussions. 



In the experiment of Strecker et al. [3] "missing soli- 
tons" are frequently observed during the time evolution 
of the train [16]. Our results strongly suggest that the 
phenomenon of "missing solitons" is related to the inter- 
mittent dynamics of individual peaks in the train under 
axial harmonic confinement. 

To simulate the finite resolution of the experimental 
detection and imaging process, we calculate the covolu- 
tion p(z) — J dz'G(z — z')p(z') of the axial density profile 
p{z) with a Gaussian G(z) having the width of a single 
soliton. As shown in Figure 5, for x = 1 the number N p 
of peaks of the smoothed density p(z) oscillates in time 
around a mean value smaller than N s . In fact, many 
peaks of the train aro too close to be seen by using a 
finite resolution of the density: that becomes dramatic 
at the turning points. Thus, the effect of finite resolu- 
tion could explain the disagreement between the exper- 
imentally observed number of solitons and its analytical 
estimate based on modulational instability. 

In Figure 5 it is also shown the behavior of N p with 
initially 10 4 atoms. Remarkably, in this case the largest 
values of N p are obtained when the train is in the middle 
of the trap: contrary to the case with initially 10 5 atoms, 
the solitons bunch in the middle and spread out at the 
turning points. Such a phenomenon can be qualitatively 
explained observing that the energy of the repulsive inter- 
action between solitons decreases by reducing the number 
of atoms [5,14]. It follows that below a critical threshold 
in the number of atoms the repulsive interaction between 
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